******************************************************************************
*										
*                              Boese & Eberhardt							
*            
*                             Facets of Democracy
*													
******************************************************************************
*                   PCDID - heterogeneous model (1959-2018)
*			    (Specification with exports/trade & pop growth)
* 	                           Low-level indices
*                Varying the cutoff: Alpha and Chi2 tests
******************************************************************************
*                       Using standardized thresholds 
******************************************************************************
*                           Created: 5th June 2024					
*                      Last Changed: 5th June 2024					
******************************************************************************
*                              Markus Eberhardt						
******************************************************************************
*                 School of Economics, University of Nottingham,			
*                  University Park, Nottingham NG7 2RD, England			
*                     email: markus.eberhardt@nottingham.ac.uk				
*                   web: https://sites.google.com/site/medevecon/			
******************************************************************************

clear all
set matsize 11000

***
*** Paths
***

global path "D:/Dropbox/Facets of Democracy"
global path "/Users/lezme/Dropbox/Facets of Democracy"
global rawpath "/Users/lezme/Dropbox"

global data "$path/Data"
global dofolder "$path/Do_files"
global otherdata "$path/Data"
global outputfolder "$path/Output"
global texfolder "$rawpath/Apps/Overleaf/Facets of Democracy/images"

global xtmgpath = "D:/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"
global xtmgpath = "/Users/lezme/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"

****
**** Get merged dataset and prepare variables
****

use "$data/madd_ert_vdem_2021.dta", clear
tsset nwbcode year

* Maddison: per capita GDP (in logs)*100 and population growth rate
gen y = ln(gdppc)*100
gen lpop = ln(pop)
gen dlpop = 100*d.lpop

* DOTS data: export/trade
gen ex = ex_trade*100

xtreg y dlpop ex v2x_libdem if year>1958, fe
* n=7,643, N=157, avg T=48.7 (min 12, max 60); 1959-2018

drop csum 
gen c=1 if e(sample)
by nwbcode: egen csum=sum(c) if c==1
gen sample= (e(sample))

sum v2x_freexp_altinf v2x_frassoc_thick v2xel_frefair v2xcl_rol v2x_jucon v2xlg_legcon if sample==1


* V-Dem: democracy dummies from high and mid-level indicators
local z "_freexp_altinf _frassoc_thick el_frefair cl_rol _jucon lg_legcon"
foreach k of local z{
	sum v2x`k' if sample==1
	scalar mean`k' = r(mean) 
	scalar sd`k' = r(sd) 
	gen st`k' = (v2x`k'-mean`k')/sd`k' if sample==1
}


* Positive values
local z "0 25 50"
foreach k of local z{
	local norm_`k' = `k'/100
	gen dem_exp_`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold Threshold > mean + .`k' sd"
	gen dem_rolaw_`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legon_`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legon_`k' "Leg. Constraints Threshold > mean + .`k' sd"
}

* Positive values
local z "125"
foreach k of local z{
	local norm_`k' = `k'/1000
	gen dem_exp_`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold > mean + .`k' sd"
	gen dem_rolaw_`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legon_`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legon_`k' "Leg. Constraints Threshold > mean + .`k' sd"
}

* Negative values
local z "50 25"
foreach k of local z{
	local norm_neg`k' = -`k'/100
	gen dem_exp_neg`k' = (st_freexp_altinf>=`norm_neg`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_neg`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_neg`k' = (st_frassoc_thick>=`norm_neg`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_neg`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_neg`k' = (stel_frefair>=`norm_neg`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_neg`k' "Clean Election Threshold Threshold > mean + .`k' sd"
	gen dem_rolaw_neg`k' = (stcl_rol>=`norm_neg`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_neg`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_neg`k' = (st_jucon>=`norm_neg`k'') & !missing(v2x_jucon)
	label var dem_jucon_neg`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legon_neg`k' = (stlg_legcon>=`norm_neg`k'') & !missing(v2xlg_legcon) 
	label var dem_legon_neg`k' "Leg. Constraints Threshold > mean + .`k' sd"
}


* Negative values
local z "125"
foreach k of local z{
	local norm_neg`k' = -`k'/1000
	gen dem_exp_neg`k' = (st_freexp_altinf>=`norm_neg`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_neg`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_neg`k' = (st_frassoc_thick>=`norm_neg`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_neg`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_neg`k' = (stel_frefair>=`norm_neg`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_neg`k' "Clean Election Threshold Threshold > mean + .`k' sd"
	gen dem_rolaw_neg`k' = (stcl_rol>=`norm_neg`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_neg`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_neg`k' = (st_jucon>=`norm_neg`k'') & !missing(v2x_jucon)
	label var dem_jucon_neg`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legon_neg`k' = (stlg_legcon>=`norm_neg`k'') & !missing(v2xlg_legcon) 
	label var dem_legon_neg`k' "Leg. Constraints Threshold > mean + .`k' sd"
}


****
**** PCDID 
**** 

estimates clear 

local i = 1

**** Select the variables, sample period, factor-augmentation 
global xvars1 		"ex dlpop"	// X variables used to construct factors 
global xvars2 		"ex dlpop"	// Xs included in the PCDID (leave blank for 'plain vanilla' model)
global factors 		"6"			// max # of factors 
global uhat 		"uhatMG"	// uhatMG: residuals from heterog regression, uhatP: pooled FE

* Need to use 2 sets as esttab can only deal with 150 models max (one commented out!!!)
* Polyarchy components
*global demo_list	"exp_neg25 exp_neg125 exp_0 exp_125 exp_25  assoc_neg25 assoc_neg125 assoc_0  assoc_125  assoc_25 fairel_neg25 fairel_neg125 fairel_0 fairel_125 fairel_25" 
*global demo_list2	"dem_exp_neg25 dem_exp_neg125 dem_exp_0 dem_exp_125 dem_exp_25 dem_assoc_neg25 dem_assoc_neg125 dem_assoc_0  dem_assoc_125  dem_assoc_25 dem_fairel_neg25 dem_fairel_neg125 dem_fairel_0 dem_fairel_125 dem_fairel_25" 

* Liberal component components
global demo_list 	"rolaw_neg25 rolaw_neg125 rolaw_0 rolaw_125 rolaw_25 jucon_neg25 jucon_neg125 jucon_0 jucon_125 jucon_25 legon_neg25 legon_neg125 legon_0 legon_125 legon_25"
global demo_list2	"dem_rolaw_neg25 dem_rolaw_neg125 dem_rolaw_0 dem_rolaw_125 dem_rolaw_25 dem_jucon_neg50 dem_jucon_neg25 dem_jucon_neg125 dem_jucon_0 dem_jucon_125 dem_jucon_25 dem_legon_neg25 dem_legon_neg125 dem_legon_0 dem_legon_125 dem_legon_25"

* Full list (for merger and graphs)
global demo_list_full	"exp_neg25 exp_neg125 exp_0 exp_125 exp_25 assoc_neg25 assoc_neg125 assoc_0  assoc_125  assoc_25 fairel_neg25 fairel_neg125 fairel_0 fairel_125 fairel_25 rolaw_neg25 rolaw_neg125 rolaw_0 rolaw_125 rolaw_25 jucon_neg25 jucon_neg125 jucon_0 jucon_125 jucon_25 legon_neg25 legon_neg125 legon_0 legon_125 legon_25"

run "$xtmgpath/xtmg.ado"

* Start of loop for different end years 
local zzz "2018"
foreach kkk of local zzz{
	global endyear 		"`kkk'"			// 2018

	* Start of loop for different start years 
	local zz "1959"
	foreach kk of local zz{
		global startyear 	"`kk'"		// 1959 

		* Start of loop for different democracy indicators 
		local z "$demo_list"
		foreach k of local z{
			global democracy 	"`k'"		// pick democracy indicator (V-Dem indices)

			* Start of factor loop: how many factors to be included in the PCDID
			forvalues f=1/$factors{
				preserve
				quietly{
					tsset nwbcode year
					keep if year>=$startyear & year<=$endyear
					gen dem_dum = dem_$democracy
					sort nwbcode year
					* Cross the threshold
					gen count_$democracy = 1 if ((dem_dum==0 & l.dem_dum==1) | (dem_dum==1 & l.dem_dum==0))
					by nwbcode: egen dem_dum_count_$democracy = sum(count_$democracy)
					replace count_$democracy = 0 if (dem_dum==0 & l.dem_dum==1)
					* Cross the threshold from below
					by nwbcode: egen trudem_dum_count_$democracy = sum(count_$democracy)
					by nwbcode: egen dem_dum_exposure_$democracy = sum(dem_dum)
					by nwbcode: egen sample_obs_$democracy = sum(sample)
					gen autoc_exposure_$democracy = sample_obs_$democracy - dem_dum_exposure_$democracy
					xtlogit dem_dum y if sample==1, fe iterate(1)
					* Treat are those which experienced crossing the threshold
					gen treat = (e(sample))
					replace treat = . if sample==0
					* Below are those which are always below the threshold
					gen below = (!e(sample) & sample==1 & !missing(dem_$democracy))
					replace below = . if sample==0
					* Always above threshold (adjust)
					gen above = (below==1 & dem_$democracy ==1)
					replace above = . if sample==0
					* Always below threshold (adjust)
					replace below = . if below==1 & dem_$democracy==1
					
					*** Estimate linear heterog regression of y on x in the control group and collect residuals
					xtmg y $xvars1 if below==1 & sample==1 & above==0, res(uhatMG)
					local below = e(N_g)
					local belowobs = e(N)
					*** Estimate linear pooled regression of y on x in the control group and collect residuals
					xtreg y $xvars1 if below==1 & sample==1, fe
					predict uhatP, e
					*** Compute residual period average in control group 
					sort year 
					by year: egen uCT = mean($uhat)
					sort nwbcode year					
					*** Extract factors from the residuals - decide which version ($uhat) to use
					xtreg $uhat if sample==1, fe
					gen t_max = e(g_max)
					gen N_max = e(N_g)
					if t_max > N_max{
						regife $uhat if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
					}
					else if  N_max > t_max {
						regife $uhat if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
					}
					drop t_max N_max
					
					*** Create factors for all countries
					forvalues l=1/`f'{
						sort year
						by year: egen cfactor`l'N = mean(cfactor`l')		
						drop cfactor`l' 
					}
					drop floading*
					
					*** Test of controls ***
					xtlogit dem_$democracy y $xvars2 ///
						if sample==1 & treat==1, fe iterate(1)
					xtmg dem_$democracy $xvars2 cfactor*N ///
						if e(sample)
					testparm $xvars2
					local chi2_stat = r(chi2)
					local chi2_p = r(p)

					*** Alpha test auxiliary regressions ***
					xtmg y uCT ///
						dem_$democracy ///
						$xvars2 ///
						if sample==1 & treat==1
					mat alphas = e(betas)
					svmat alphas, name(alphas)
					gen alpha01 = alphas1-1
					replace alpha01 = . if (alphas2)==0
					mkmat alpha01, mat(alpha01)
					set seed 280115
					reg alpha01
					mat a1 = e(b)
					mat vv1 = e(V)
					scalar alphaest1 = a1[1,1]
					local alpha_est1 = alphaest1
					scalar vv1s = vv1[1,1]
					scalar se1s = sqrt(vv1s)
					scalar alphat1 = alphaest1/se1s
					local alpha_tstat1 = alphat1
					
					*** Run PCDID model ***
					eststo: xtmg y dem_$democracy $xvars2 ///
						cfactor*N ///
						if sample==1 & treat==1, res(resid) robust
					mat betas = e(betas)
					estadd scalar treated = e(N_g)
					estadd scalar obs = e(N)
					local treated = e(N_g)
					local obs = e(N)
					mat b1 = e(b)
					mat v1 = e(V)
					scalar b_1 = b1[1,1]
					scalar se_1 = sqrt(v1[1,1])	
					scalar mse = e(sigma)
					estadd scalar mse = e(sigma)
					capture xtcd resid if e(sample)
					capture scalar cdtest = r(pesaran)
					capture estadd scalar cdtest = r(pesaran)
					estadd scalar fac = `f'
					estadd scalar x_test_chi2 = `chi2_stat'
					estadd scalar x_test_p = `chi2_p'
					estadd scalar alpha_b1 = `alpha_est1'
					estadd scalar alpha_t1 = `alpha_tstat1'
					estadd scalar controlN = `below'
					estadd scalar controlobs = `belowobs'
					estadd scalar final_year = `kkk'
					estadd scalar start_year = `kk'
					est store model`i'
					
					keep if sample==1 & treat==1
					tabstat nwbcode if sample==1 & treat==1, by(wbcode) save
					tabstatmat nwbcodestat, nototal
					tabstat year if sample==1 & treat==1, by(wbcode) stats(min max) save
					tabstatmat yearstat, nototal
					tabstat dem_dum_count_$democracy trudem_dum_count_$democracy dem_dum_exposure_$democracy ///
						sample_obs_$democracy autoc_exposure_$democracy, by(wbcode) stats(mean) save
					tabstatmat exposure, nototal
					drop _all
					svmat nwbcodestat, name(nwbcode)
					rename nwbcode1 nwbcode
					svmat yearstat, name(year)
					rename year1 startyear
					rename year2 endyear
					svmat exposure, name(count)
					rename count1 flips
					rename count2 democratisations
					rename count3 exposure
					rename count4 total_obs
					rename count5 nonexposure
					svmat betas, name(beta)
					gen wbcode=""
					replace wbcode="AFG" if nwbcode==	2
					replace wbcode="AGO" if nwbcode==	3
					replace wbcode="ALB" if nwbcode==	5
					replace wbcode="ARE" if nwbcode==	7
					replace wbcode="ARG" if nwbcode==	8
					replace wbcode="ARM" if nwbcode==	9
					replace wbcode="AUS" if nwbcode==	12
					replace wbcode="AUT" if nwbcode==	13
					replace wbcode="AZE" if nwbcode==	14
					replace wbcode="BDI" if nwbcode==	15
					replace wbcode="BEL" if nwbcode==	17
					replace wbcode="BEN" if nwbcode==	18
					replace wbcode="BFA" if nwbcode==	19
					replace wbcode="BGD" if nwbcode==	20
					replace wbcode="BGR" if nwbcode==	21
					replace wbcode="BHR" if nwbcode==	22
					replace wbcode="BIH" if nwbcode==	24
					replace wbcode="BLR" if nwbcode==	25
					replace wbcode="BOL" if nwbcode==	29
					replace wbcode="BRA" if nwbcode==	30
					replace wbcode="BRB" if nwbcode==	31
					replace wbcode="BWA" if nwbcode==	36
					replace wbcode="CAF" if nwbcode==	37
					replace wbcode="CAN" if nwbcode==	38
					replace wbcode="CHE" if nwbcode==	39
					replace wbcode="CHL" if nwbcode==	40
					replace wbcode="CHN" if nwbcode==	41
					replace wbcode="CIV" if nwbcode==	42
					replace wbcode="CMR" if nwbcode==	43
					replace wbcode="COG" if nwbcode==	45
					replace wbcode="COL" if nwbcode==	46
					replace wbcode="COM" if nwbcode==	47
					replace wbcode="CPV" if nwbcode==	48
					replace wbcode="CRI" if nwbcode==	49
					replace wbcode="CUB" if nwbcode==	51
					replace wbcode="CYP" if nwbcode==	52
					replace wbcode="CZE" if nwbcode==	53
					replace wbcode="DEU" if nwbcode==	55
					replace wbcode="DJI" if nwbcode==	56
					replace wbcode="DNK" if nwbcode==	58
					replace wbcode="DOM" if nwbcode==	59
					replace wbcode="DZA" if nwbcode==	60
					replace wbcode="ECU" if nwbcode==	61
					replace wbcode="EGY" if nwbcode==	62
					replace wbcode="ESP" if nwbcode==	64
					replace wbcode="EST" if nwbcode==	65
					replace wbcode="ETH" if nwbcode==	66
					replace wbcode="FIN" if nwbcode==	67
					replace wbcode="FRA" if nwbcode==	69
					replace wbcode="GAB" if nwbcode==	71
					replace wbcode="GBR" if nwbcode==	72
					replace wbcode="GEO" if nwbcode==	73
					replace wbcode="GHA" if nwbcode==	74
					replace wbcode="GIN" if nwbcode==	76
					replace wbcode="GMB" if nwbcode==	77
					replace wbcode="GNB" if nwbcode==	78
					replace wbcode="GNQ" if nwbcode==	79
					replace wbcode="GRC" if nwbcode==	80
					replace wbcode="GTM" if nwbcode==	83
					replace wbcode="HKG" if nwbcode==	87
					replace wbcode="HND" if nwbcode==	89
					replace wbcode="HRV" if nwbcode==	91
					replace wbcode="HTI" if nwbcode==	92
					replace wbcode="HUN" if nwbcode==	93
					replace wbcode="IDN" if nwbcode==	95
					replace wbcode="IND" if nwbcode==	96
					replace wbcode="IRL" if nwbcode==	97
					replace wbcode="IRN" if nwbcode==	98
					replace wbcode="IRQ" if nwbcode==	99
					replace wbcode="ISL" if nwbcode==	100
					replace wbcode="ISR" if nwbcode==	101
					replace wbcode="ITA" if nwbcode==	102
					replace wbcode="JAM" if nwbcode==	103
					replace wbcode="JOR" if nwbcode==	104
					replace wbcode="JPN" if nwbcode==	105
					replace wbcode="KAZ" if nwbcode==	106
					replace wbcode="KEN" if nwbcode==	107
					replace wbcode="KGZ" if nwbcode==	108
					replace wbcode="KHM" if nwbcode==	109
					replace wbcode="KOR" if nwbcode==	112
					replace wbcode="KWT" if nwbcode==	113
					replace wbcode="LAO" if nwbcode==	115
					replace wbcode="LBN" if nwbcode==	116
					replace wbcode="LBR" if nwbcode==	117
					replace wbcode="LBY" if nwbcode==	118
					replace wbcode="LKA" if nwbcode==	120
					replace wbcode="LSO" if nwbcode==	121
					replace wbcode="LTU" if nwbcode==	122
					replace wbcode="LUX" if nwbcode==	123
					replace wbcode="LVA" if nwbcode==	124
					replace wbcode="MAR" if nwbcode==	126
					replace wbcode="MDA" if nwbcode==	128
					replace wbcode="MDG" if nwbcode==	129
					replace wbcode="MEX" if nwbcode==	132
					replace wbcode="MLI" if nwbcode==	135
					replace wbcode="MLT" if nwbcode==	136
					replace wbcode="MMR" if nwbcode==	137
					replace wbcode="MNE" if nwbcode==	138
					replace wbcode="MNG" if nwbcode==	139
					replace wbcode="MOZ" if nwbcode==	140
					replace wbcode="MRT" if nwbcode==	141
					replace wbcode="MUS" if nwbcode==	143
					replace wbcode="MWI" if nwbcode==	144
					replace wbcode="MYS" if nwbcode==	145
					replace wbcode="NAM" if nwbcode==	146
					replace wbcode="NER" if nwbcode==	148
					replace wbcode="NGA" if nwbcode==	149
					replace wbcode="NIC" if nwbcode==	150
					replace wbcode="NLD" if nwbcode==	151
					replace wbcode="NOR" if nwbcode==	152
					replace wbcode="NPL" if nwbcode==	153
					replace wbcode="NZL" if nwbcode==	156
					replace wbcode="OMN" if nwbcode==	158
					replace wbcode="PAK" if nwbcode==	159
					replace wbcode="PAN" if nwbcode==	160
					replace wbcode="PER" if nwbcode==	161
					replace wbcode="PHL" if nwbcode==	162
					replace wbcode="POL" if nwbcode==	165
					replace wbcode="PRK" if nwbcode==	168
					replace wbcode="PRT" if nwbcode==	170
					replace wbcode="PRY" if nwbcode==	171
					replace wbcode="QAT" if nwbcode==	176
					replace wbcode="RUS" if nwbcode==	179
					replace wbcode="RWA" if nwbcode==	180
					replace wbcode="SAU" if nwbcode==	181
					replace wbcode="SDN" if nwbcode==	183
					replace wbcode="SEN" if nwbcode==	184
					replace wbcode="SGP" if nwbcode==	185
					replace wbcode="SLE" if nwbcode==	187
					replace wbcode="SLV" if nwbcode==	188
					replace wbcode="STP" if nwbcode==	195
					replace wbcode="SVK" if nwbcode==	198
					replace wbcode="SVN" if nwbcode==	199
					replace wbcode="SWE" if nwbcode==	200
					replace wbcode="SWZ" if nwbcode==	201
					replace wbcode="SYC" if nwbcode==	203
					replace wbcode="SYR" if nwbcode==	204
					replace wbcode="TCD" if nwbcode==	205
					replace wbcode="TGO" if nwbcode==	206
					replace wbcode="THA" if nwbcode==	207
					replace wbcode="TJK" if nwbcode==	208
					replace wbcode="TKM" if nwbcode==	209
					replace wbcode="TTO" if nwbcode==	214
					replace wbcode="TUN" if nwbcode==	215
					replace wbcode="TUR" if nwbcode==	216
					replace wbcode="TZA" if nwbcode==	220
					replace wbcode="UGA" if nwbcode==	221
					replace wbcode="UKR" if nwbcode==	222
					replace wbcode="URY" if nwbcode==	223
					replace wbcode="USA" if nwbcode==	224
					replace wbcode="UZB" if nwbcode==	225
					replace wbcode="VEN" if nwbcode==	228
					replace wbcode="VNM" if nwbcode==	229
					replace wbcode="YEM" if nwbcode==	235
					replace wbcode="ZAF" if nwbcode==	238
					replace wbcode="ZMB" if nwbcode==	240
					replace wbcode="ZWE" if nwbcode==	241
					order nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
					keep nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
					gen regime_def = "$democracy"
					local startyearx "$startyear"
					local endyearx "$endyear" 
					local democracyx "$democracy" 
					save "$outputfolder/PCDID_Maddison2024_Drilling_`startyearx'_`endyearx'_st_multiple_`democracyx'_diagnostics_`f'.dta", replace
					
					local i = `i'+1
					
						}
					* end of quietly
					
				di in ye _n "*****************************************************************"
				di in ye "*** Single regime cutoff defined by " in gr "dem_$democracy" 
				di in ye "*****************************************************************"
				di in ye "*** PCDID (MG) estimate with" in gr " `f' " in ye "common factor(s) included      ***"
				di in ye "*****************************************************************"
*				di in ye "Bai & NG: PCp1 " %5.3f PC1_`f' 
				di in ye "*****************************************************************"
				di in ye "Nominal sample period: " in gr "$startyear" in ye " to " in gr "$endyear "
				di in ye "resulting in " in gr "`treated'" in ye " treated countries with " ///
					in gr "`obs'" in ye " observations. "
				di in ye "MSE: " in gr %5.3f mse
				di in ye "Control sample: " in gr "`below'" in ye " countries with " in gr "`belowobs'" in ye " obs." 
				di in ye "Rob. Beta: " in gr %5.3f b_1 in ye "  SE: " in gr %5.3f se_1 ///
					in ye "  t-ratio: " in gr %5.3f b_1/se_1 in ye "  p: " in gr %5.3f 2*(normal(-abs(b_1/se_1)))
				di in ye "*****************************************************************"
				di in ye "Alpha Test: " in gr  %5.3f alphaest1 in ye " t-stat: " in gr  %5.3f alphat1
				di in ye "*****************************************************************"
				di in ye "*** Test of controls " in gr %5.2f `chi2_stat' in ye " (" in gr %5.3f `chi2_p' in ye ").      ***"
				di in ye "*****************************************************************"
				local ii = `i'-1
				di in red "`ii'"
				
				restore	
				
			}
			* end of the factor loop
		}
		* end of dem indicator loop
	}
	* end of start year loop
}
* end of end year loop

local mytitle "PCDID"
local esttabformat  "nodepvars lines compress label nogaps numbers star(* 0.10 ** 0.05 *** 0.01)"
local esttabstats "stats(start_year final_year fac treated obs controlN controlobs x_test_chi2 x_test_p alpha_b1 alpha_t1, fmt(0 0 0 0 0 0 0 2 2 3 3))"
local esttabkeep "keep($demo_list2 $xvars2) order($demo_list2 $xvars2)"
esttab model* using "$outputfolder/20240605_PCDID_1959_2018_Maddison_Drilling_Low_level_Lib_multiple_Diagnostics.csv", b(3) se(3) abs br ///
	`esttabformat' `esttabkeep' `esttabstats' style(tab) replace nodep

** Change the name according to which sub-components were studied
*20240605_PCDID_1959_2018_Maddison_Drilling_Low_level_Poly_multiple_Diagnostics
*20240605_PCDID_1959_2018_Maddison_Drilling_Low_level_Lib_multiple_Diagnostics
	
exit	
